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Abstract 

We demonstrate that stabilization of solitons of the multidimensional Schrodinger 
equation with a cubic nonlinearity may be achieved by a suitable periodic control 
of the nonlinear term. The effect of this control is to stabilize the unstable solitary 
waves which belong to the frontier between expanding and collapsing solutions and 
to provide an oscillating solitonic structure, some sort of breather-type solution. 
We obtain precise conditions on the control parameters to achieve the stabiliza- 
tion and compare our results with accurate numerical simulations of the nonlinear 
Schrodinger equation. Because of the application of these ideas to matter waves 
these solutions are some sort of matter breathers. 
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waves 
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1 Introduction 



The nonlinear Schrodinger equation (NLSE) in its many versions is one of the 
most important models of mathematical physics, with applications to different 
fields such as plasma physics, nonlinear optics, water waves and biomolecular 
dynamics, to cite only a few cases. In many of those examples the equation 
appears as an asymptotic limit for a slowly varying dispersive wave envelope 
propagating in a nonlinear medium [1,2]. 
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A new burst of interest on problems modelized by nonlinear Schrodinger equa- 
tions has been triggered by the achievement of Bose-Einstein condensation 
using ultracold neutral bosonic gases [3,4]. 

In particular in Bose-Einstein condensation one of the key ingredients of the 
achievement of condensation with alkali gases is the trapping of the neutral 
atoms which are cooled down below the transition temperatures. Although 
different trapping techniques have been used in practice, the most commonly 
used are the magnetic traps which are modelized by external confining forces 
acting on the system. 

There are specific types of Bose-Einstein condensates such as those made 
of Lithium [5,6] in which the interactions between the atoms are attractive. 
Mathematically this implies that for spatial dimensionalities n = 2, 3 collaps- 
ing solutions are possible as it is well known in the framework of the usual 
studies of nonlinear Schrodinger equations [2]. Thus, although in principle one 
could think that nonlinearity and dispersion could balance mutually and a self- 
sustained solitonic structure might exist this is only true for n = 1. In fact, 
it was soon found experimentally [5,6] and theoretically supported that such 
a situation is unstable and leads to collapse [2,7]. Later studies, taking into 
account more elaborate models than the simpler NLS mean field models, led 
to the understanding that the occurrence of collapse during the condensation 
process would limit the size of an attractive condensate [9,10]. 

Thus, the only confirmed way to observe solitonic states in Bose-Einstein 
condensates (matter-wave solitons) with negative scattering length involves 
the elimination of the trap in one direction as proposed in [11,12] and found 
experimentally in [13,14]. 

However, the possibility of using Feschbach resonances to control the scatter- 
ing length [15,16,17] has provided a way to study large negative scattering 
length condensates and collapse processes in detail [18,19]. This control al- 
lows, with certain experimental limitations, to modulate appropriately the 
nonlinear term. 

But this possibility opens new ways for the generation and observation of dif- 
ferent types of matter-wave solitons not yet fully explored. One of the most 
striking ones is the possibility of generating trapless trapped Bose-Einstein con- 
densate solitons. The idea is that (oscillating) bound states might be obtained 
by combining cicles of positive and negative scattering length values so that 
after an expansion and contraction regime the condensate would come back to 
the initial state. In this way some sort of pulsating trapped condensate, i.e., 
a breather, would be obtained. 

A rigorous study of this intuitive idea is necessary in order to make precise 
predictions and find which are the precise parameter values to be used to ob- 
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tain such breathers. This is the purpouse of this paper in a general framework 
of nonlinear Schrodinger systems. This idea has been explored in two previous 
papers [20,21] but here we improve the understanding of the phenomenon and 
correct mistakes contained in the previous analysis. 

Another field of applications of these ideas is nonlinear Optics, there a full 
control of the nonlinear term is not possible but, because of technical limita- 
tions, only piecewise constant values for the nonlinear coefficient can be easily 
generated experimentally. 

/,From the mathematical point of view what we would like is to stabilize solu- 
tions close to the stationary unstable solitons of the cubic NLSE by choosing 
appropriate controls of the nonlinear term. To our knowledge this problem has 
not been considered previously in the mathematical literature. 

This paper is organized as follows: First, in Sec. 2 we present the model equa- 
tions as they appear in one specific application of the model. Next, in Sec. 3 
we reduce the nonlinear model to a system of ordinary differential equations 
by using the so-called moment method. In Sec. 4 we analyze two-dimensional 
systems and compare the analytical predictions of the moment method with 
numerical simulations of the full partial differential equations. We obtain pre- 
cise conditions for the existence of breathers and discuss the importance of the 
initial data. In Sec. 5 we analyze the three-dimensional case and explain why 
moment (or related) equations fail to describe the dynamics of the spherically 
symmetric case. We consider the case of systems in specific external potentials 
and give some stabilization conditions. Finally, in Sec. 6 we summarize our 
conclusions and compare our results with the previous findings of Refs. [20,21]. 



2 The model 



In this work we study systems ruled by the NLSE with a cubic nonlinearity. In 
the framework of recent problems in Bose-Einstein condensation this equation 
appears as the model for the mean field dynamics of a boson system described 
by a single wavefunction \1> in the zero-temperature limit. The resulting NLSE 
equation is called sometimes the Gross-Pitaevskii equation (GPE) which is 
[22] 
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where A = X)f=i d 2 / dr 2 is the three-dimensional Laplacian, V(r) — ^ ■ . , r 

is the external potential (the so called "trap") which confines the condensate 
and a is the s-wave scattering length for the binary collisions within the con- 
densate. 
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It is more convenient to work with a new set of adimensional quantities defined 
r-j/ao, t = t/T, ip(x., t) = ^(r, t)\Jo,q/N, where a = \j h/mu, T — 1/uj 



as x 



and N — J \ty\ 2 d 3 r is the number of particles in the condensate. Then Eq. (1) 
reads as the next NLSE 
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Z 3=1 



1>, 



(2) 



where g = AitN a /do and / |-^| 2 (i 3 x = 1. In this paper we consider a system 
evolving without external potential along one or more dimensions. In Bose- 
Einstein condensation this corresponds to a condensate generated in a trap 
and then the trapping is eliminated. This leads to different particular cases of 
Eq. (2). When the potential is removed in all directions we obtain 



.dip 



ip. 



(3) 



This equation holds for three-dimensional and quasi two-dimensional systems 
as will be discussed in detail in Sec. 5. Eq. (3) has the typical form of the 
NLSE with cubic power nonlinearity [2] and has been extensively studied. In 
this work we will analyze the possibilities of controlling the behavior of the 
stationary solutions of Eq. (3) by letting g be a time-dependent real function. 
Similar possibilities for nonconstant g (in that context z <-> t) but restricted 
to n = 2 arise in the context of coherent light propagation in nonlinear Kerr 
media in the paraxial approximation. 



3 Moment method 



In this section we build a theory which will allow us to reduce the dynamics 
from the complex model given by Eq. (3) to a simpler one whose analysis will 
provide insight on the basic features of the phenomenon. 

Let us first consider radially symmetric solutions of Eq. (3) u = u(r,t), r = 
(E"=i^) 1/2 , satisfying 

^=-^#(H^W*)i«i 2 «- (4) 

dt 2r n ~ l dr\ dr ) y ' 

where n is the spatial dimensionality, in our case n — 2, 3. To get information 
on the solutions of Eq. (4) we will use the moment method [23,24,25,26,27]. 
This method proceeds by analyzing the evolution of several integral quantities 
[23,24,25] related with Eq. (4). 
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1^) = J \u\ 2 d n x (5a) 
I 2 (t) = J \u\ 2 r 2 d n x, (5b) 

h(t) = i [(u^-u*^)rd n x, (5c) 




/ 4 (t) = -/ (\Vu\ 2 + ^g(t)\u\^d n x, (5d) 
I 5 (t) = jj \u\ 4 d n x, (5e) 

where cPa; = r n_1 <ir<if2. With our scaling for *0, the first one satisfies h = 1, for 
all t. The remaining ones are related physically to the width, radial momentum 
and energy of the wave packet. In what follows we assume that the initial data 
are such that all Ij are initially well-defined [28]. Some algebra leads to 



£='•• « & » 

f=4/„ (6b, 

dh _ n-2dl 5 ^ dg j 
dt n dt dt 5 ' 

dh _ mr2- r d\u\ 4 d^gu x _ _nhh , 

All the momenta follow closed evolution laws except for In specific 

situations the moment equations provide completely closed equations for the 
evolution of all the Ij [25,27,28] but this is not our case. Extending the number 
of momenta which are included in the calculation does not help to obtain a 
closed set of equations. 

Thus, to close the system we will take 

argw = I 3 r 2 /AI 2 , (7) 

as used in [26]. Physically, this corresponds to approximating the phase of u by 
the spherical wavefront which best fits the distribution. A rigorous justification 
of this choice is possible [8] for the case n = 2 and when the initial data is the 
so-called Townes soliton to be presented later. 

When solutions with phase given by Eq. (7) are considered Eqs. (6) are closed 
and have several (positive) invariants under time evolution given by 



Q 1 = 2(I 4 -gI 5 )I 2 -I 2 3 /4, (8a) 
Q 2 = 2I^ /2 I 5 . (8b) 
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With the help of these quantities, Eqs. (6) can be reduced to a single equation 
for hit), which is 



If we were able to solve Eq. (9) then the use of Eqs. (6) would allow us to track 




the evolution of the other ones. To simplify Eq. (9) we define X{t) = ^^{t), 
which is the wave packet width, and substituting it into (9) we get 



4 Two-dimensional systems 



4-1 Rigorous analysis of the moment equations 



In physical situations both problems with n = 2 and n = 3 arise. The situ- 
ation when a Bose-Einstein condensate is tightly confined along a particular 
direction leads to a quasi two-dimensional system such as the ones obtained in 
Ref. [29]. From now on we consider the two-dimensional model to be valid in 
that situation although a more detailed analysis will be made in Sec. 6. Thus, 
when n = 2 we obtain 

X ~ ^3 =xs- ( n ) 

This equation appears also in relation with the classical motion of an atom 
near an infinite straight wire with oscillating charge (the Paul trap) [30,31,32] 
and also as an approximate model for beam propagation in layered Kerr media 
[33]. 

Let us consider the case where p is continuous and T-periodic and look for 
positive periodic solutions and bound states, that is, bounded (nonperiodic) 
positive solutions without collapse. If p(t) > all solutions escape to infinity 
whereas if pit) < all solutions collapse. In fact, theorem 2.1 of [31] implies 
that a necessary condition for existence of bound states is 

p=U T p<0. (12) 



T Jo 



Let pit) be parametrized as p(t) = a + (5a{t) with a = 0. We can fix a M = 
maxa(t) = 1 without loss of generality. A direct consequence of the latter 



6 



result is that if a > bound states cannot exist, therefore a must be negative. 
Also, p(t) must change sign, so for existence of bound states it is necessary 
that 

a + \P\>0. (13) 



Let us introduce the positive small parameter e = — 2a/3\P\. We will now 
prove that there exists e > (depending only on a) such that Eq. (11) has a 
Lyapunov- stable periodic solution if < e < e . Moreover, there is an infinite 
number of quasiperiodic solutions and a sequence of subharmonics with min- 
imal periods tending to +oo. To prove this affirmation let us note that when 
a(t) = cos fit, Theorem 1 in [32] asserts that if e is small enough there is a pe- 
riodic solution of twist type [34]. However, the proof is still valid for a general 
continuous and T-periodic function with a = 0. Moser Twist Theorem [35] im- 
plies that a solution of twist type is Lyapunov-stable as a consequence of the 
existence of invariant curves (quasiperiodic solutions) around it. Finally, the 
existence of subharmonics of any order is a consequence of Poincare-Birkhoff 
Theorem. 

This result has interesting physical consequences. The existence of periodic so- 
lutions is directly related to the existence of pulsating breathers. The bounded 
solutions of the theorem would correspond to solutions with quasiperiodic or 
chaotic oscillations in the width. In fact, Mather sets and Smale's horseshoes 
appear in the Poincare map. The latter correspond to chaotic oscillations 
whereas Mather sets are Cantorian structures that are more complicated than 
the quasiperiodic ones. 

We have just seen that a necessary condition for the existence of bounded 
solutions is that a+ \/3\ > with a < 0. Moreover the previous result indicates 
that the amplitude of the oscillating term must be of the same order as the 
amplitude of the non-oscillating term. 

To find particular types of these breathers let us take g(t) as g(t) = go + 
gi cos(fii) so that p(t) = a + (3 cos(fit) with a = Q\ + g^Q2 and f3 = g\Q2- 
An useful tool to get periodic solutions is the stability diagram of Eq. (11). 
We look for the region of parameters for which there exist stable periodic 
solutions. To do this we reduce the number of parameters by defining a new 
time as i = Qt. Then Eq. (11) is 

X" = a + ^ C 3 QS(f) , (14) 

where the prime denotes derivative with respect to t, a — a/Q 2 and f3 = 
(3/n 2 . Note that we use the word stable with two different meanings: the more 
mathematical one to indicate Lyapunov-stability and the more physical one to 
indicate existence of periodic solutions of the nonlinear Schrodinger equation. 
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Our objective is to find the region of parameters (a,/3) in which there exist 
initial conditions X = X(t = 0), X' Q = X'(t = 0) = for which Eq. (14) has 
a stable periodic solution. We already know that our search can be restricted 
to the region given by the conditions a < and a > —\/3\. To calculate the 
stabilization region we first solve numerically Eq. (14). Then we take into 
account the results that finding a periodic solution of Eq. (14) is equivalent to 
finding a solution which verifies the Neumann conditions X'(0) = = X'(2ir). 
The reason is that by extending in an even way this solution, we obtain a 
47r-periodic solution that, of course, is bounded and periodic and this implies 
that there is a 27r-periodic solution [36]. To obtain solutions that verify the 
Neumann conditions, given a pair of parameters (a, (3), we look for two initial 
conditions X^o and X2 t o for which X[ {2-k)X' 2 {2'k) < 0. Then, by continuity 
with respect to the initial conditions it is known that there exists one initial 
condition X , between Xi$ and X 2i0 , for which the solution of the equation 
satisfies X'(2tt) = 0. We will use two more observations. First of them is 
that if X{t) is a solution of the equation with parameters (a,/3), then Y(t) = 
cX(t) with c any arbitrary positive constant is also a solution with parameters 
(<5c 4 ,/3c 4 ). Therefore, by moving c, all the points of the line y = (f3/a)x are 
obtained and this is nothing but the line which passes through the points 
(at,P) and (0,0). So, to scan the region of parameters we can, for example, 
to pay attention to an specific value of a and then change continuosly the 
parameter (3. The second observation is that if X(t) is a periodic solution 
corresponding to a specific choice of (a,f3) then Y(t) = X(t + tc) is also a 
solution for parameters (a, -(3) since cos(t + n) = — cos(f). Therefore, the 
stability region is symmetric and we can restrict the search to (3 > 0. With all 
these considerations the results obtained show that the region of parameters 
for which there exist initial conditions leading to stable periodic solutions is 
given by Eq. (13) so this condition seems to be not only necessary but also a 
sufficient one. In Fig. 1 this stability region is plotted. 



4-2 Comparison with full numerical simulations of the NLSE 

In this subsection we will compare the predictions based on the ordinary dif- 
ferential equation (11) with the simulations of the full nonlinear Schrodinger 
equation (3). 

First of all we have to choose appropriate initial data uq = u(r,t = 0) in 
order to solve numerically the NLSE. Since we are going to try to stabilize the 
solutions it is reasonable to start trying to stabilize the so-called stationary 
solutions. So, let us consider stationary solutions of Eq. (3) when g = g s is 
constant, given by w(x, t) = e l/i '<3> M (x). Here $ M (x) verifies 

A$ M -2^-2(? s |<I> M | 2 <I> M = 0. (15) 
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Fig. 1. Stability diagram for Eq. (14) as a function of the parameters (a, ft). 



As it is precisely stated in [2], when g s is negative, for each positive /x there 
exists only one solution of Eq. (1) which is real, positive and radially symmet- 
ric and for which (/ \^^\ 2 d a x) 1 ^ 2 has the minimum value between all of the 
possible solutions of Eq. (15). Moreover, the positivity of // ensures that this 
solution decays exponentially at infinity. This solution is called the ground 
state and, in two dimensions, is known as the Townes soliton. We will denote 
this solution as R^ij-) and satisfies the equation 

AR, - 2/ii^ - 2g s Rl = (16) 

and the boundary conditions 

r lim^(r) = 0, #„(0) = 0. (17) 

Fixed a value of //, the norm and width of i? M are given by i]^ = (J \R^\ 2 d 2 x) 1 / 2 
and X M = (/ \R IJ \ 2 r 2 <Px) 1 l' 1 respectively. Applying scaling trasformations to 
Rn(r) it is possible to build stationary solutions having the same shape and 
norm as the Townes soliton but different widths 

R,(r) = v 1/2 Ri(» 1/2 r), X, = X,/^ 2 . (18) 
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The equation which verifies the normalized soliton R^N^r) = _R M (r)/r] M is 

AR„ tN - 2^ N - 2g s r,lRl N = 0. (19) 

We have taken as initial condition u(r,t = 0) = R^n{t) for some fi (therefore 
the initial width is fixed) and g s values. To obtain the shape R^{r) we have 
used a shooting method [37] to solve Eq. (16). The idea is to rewrite Eq. (16) 
as a dynamical system and impose that the solution of such a system also 
verifies the boundary conditions (17). In order to avoid the singularity which 
appears at the point r = we have solved the indetermination by means of a 
series expansion of R^ around r = 0. 

/From the theory of nonlinear Schrodinger equations it is known that the 
Townes soliton is unstable, i.e. small perturbations of this solution lead to 
either expansion of the initial data or blow-up in finite time. Thus, following 
the analogy with the stabilization of unstable fixed points in finite dimensional 
dynamical systems by fast variation of a parameter, we will try to stabilize 
this unstable but anyway stationary solution of Eq. (3). 

Using the predictions of Eq. (11), in order to get stable periodic solutions, we 
have to choose go such that a = Q± + goQ2 is negative, that is, the following 
relation must hold 

90 < (20) 

When only such a go-term is present the solutions collapse as predicted by Eq. 
(11). So, the gx-term is necessary in order to arrest the collapse and achieve 
the stabilization of the soliton given by u . 

Thus, we have first taken as initial data the Townes soliton obtained for g s = 
—0.5 and ji = 0.5 whose width is X = 1.09 and its norm satisfies g s vfi=o .5 = 
—5.85. This value is precisely the critical nonlinearity for collapse below which 
the collapse occurs. Let us note that only for the Townes soliton it is satisfied 
that the collapse threshold in Eq. (3), when g is a constant, is just the value 
of the nonlinearity in Eq. (19) [2,8]. 

In all the simulations of Eq. (3) we have made there appear two different types 
of behaviors. First, there is a small-amplitude fast oscillation due to the term 
cos(fit) which is present in both the full model Eq. (3) and the reduced ODE 
system. Secondly, there is a "slow" dynamics due to the internal dynamics 
of the system which sometimes does not appear in the ODE model. When 
the results from the ODE (11) and the full PDE are compared qualitatively 
similar dynamics are observed in many situations. For instance, in Fig. 2 we 
plot the results of the simulation for parameter values go = —2-k, gi = 8ir and 
Q = 40. To solve numerically Eq. (3) we have used a pseudo-spectral method 
combined with a second order split-step method to advance in time [38,39,40]. 
We have also compared this simulation with the one obtained with a fourth 
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Fig. 2. Results of numerical simulations of Eq. (3) showing stabilization of the initial 
data u(x, 0) = i? Mi jv(x) with ^ = 0.5 and g s = —0.5 (Xq = 1.09) for parameter 
values go = —2ir, g\ = 8ir, = 40. (a-c) Plots of |n(x, t)| 2 on the spatial region 
[—2,2] x [—2,2] corresponding to times of maximum (a,c) and minimum (b) width 
of the solution, (d) Evolution of the width X{t) and of the amplitude max \u(-, t)\. 
(e) Evolution of the norm i]^(t). 

order split-step method and found both results in full agreement. 

It is important to indicate that it is absolutely necessary to incorporate an 
absorbing potential at the frontier of the simulation region to avoid possible 
interferences between the fraction of the wave packet moving outside (i.e. a 
possible non-trapped fraction of initial data) and the fraction of the initial 
data which remains trapped. Typical grid sizes were 256 x 256 and the time 
step was dt = 0.001. All results were tested on different grid sizes and changing 
time steps. 
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Fig. 3. Results of numerical simulations of Eq. (3) and Eq. (11) showing stabilization 
of the initial data u(x, 0) = i?u,iv( x ) with jjl = 0.5 (Xq = 1.09) and g s = —0.5 for 
parameter values go = —2ir, g± = 8tt, f2 = 40. The evolution of the width X(t) 
according to Eq. (3) (solid) and to Eq. (11) (dashed) is shown. 

Clearly, the system is trapped with a fast modulation in the nonlinearity {g\- 
term) and the collapse process is inhibited. Also, there is a fast oscillation 
with the same frequency Q as the one of the oscillating term, (beyond the 
resolution of the plot), and a slow oscillation due to the proper dynamics of 
the system. 

In Fig. 2(e) it is observed that the norm of the solution decreases in time. This 
is caused by the action of the absorbing potential on the outgoing way. 

The next step in our study is to compare this stabilization in the NLSE with 
the differential equation for the width. In Fig. 3 we see that the quantitative 
predictions of Eq. (11) are not very precise as what concerns the amplitude 
and frequency of the slow oscillation. In any case, at least the simple ODE 
model predicts correctly the trapping of the solution for these parameters. 
These results imply that the predictions of the ODE system must be taken 
only as qualitative indications of the possible dynamics. 

Also, there exist parameter (g , g 1 and fi) choices which lead to stabilization 
in Eq. (11) but not in Eq. (3). For example, according to Eq. (11) and the 
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Fig. 4. Results of numerical simulations of Eq. (11) showing stabilization of the 
initial data w(x, 0) = i? Mi jv(x) with u = 0.5 and g s = —0.5 (X$ = 1.09) for parameter 
values go = —6.4, gi = 33.4, Q, = 20. The evolution of the width X(t) is plotted. 

diagram of Fig. 1 the parameters g = —6.4, g 1 = 33.4 and VI = 20 stabilize 
the initial soliton as shown in Fig. 4. However, this stabilization does not 
happen when we solve numerically Eq. (3), instead, after a few oscillations 
their dynamics are drastically different. 



4-3 Stabilization of other initial data 

Although a Townes soliton is a natural object to stabilize in the framework 
of Eq. (3) we may wonder (i) if the procedure described above does stabilize 
other type of initial data and (ii) if the stabilization is appropriately described 
byEqs. (11). 

From the point of view of applications of Bose-Einstein condensation the pos- 
sibility of stabilizing other initial data is important since it is not clear how 
a Townes soliton could be generated in real experiments. On the other hand, 
other initial data such as Thomas-Fermi type solutions or gaussians are much 
more natural and easier to obtain. This, together with the fact that the usual 
time dependent variational method [12] is usually developed for gaussian pro- 
file functions has lead to some interest on the possibility of trapping gaussian 
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initial data (which is the case studied in Refs. [20,21]). 



Therefore, let us take u as 



e - 2 /2X 2 

«o(r) = 7 = y : ' ( 21 ^ 



for which 



v = (^j \u \ 2 d 2 x^ = 1, (22a) 
\ 1/2 

M |V 2 rf 2 a;J . (22b) 



In this case we get for the invariants the values Q± — 1 and Q 2 = 1/2% 
independently of the width X . So, the equation for the width is 

~ l + g(t)/2% _ p(t) 

X ~ = xs- ^ 

which corresponds to the evolution equation of Refs. [20,21]. When g is a 
constant this equation predicts a threshold value for collapse of g — —2% 
which is, indeed, below the real value which corresponds to the Townes soliton 
(9sVfi = —5.85). This means that when a Gaussian function is taken as initial 
data, Eq. (23) does not describe accurately the region of trapping, because for 
— 2% < g < —5.85 this equation predicts expansion of the initial data while 
the real dynamics of the partial differential equation is collapsing. 

In Fig. 5 the evolution of a Gaussian with Xo = 1.09 and parameter values 
g = —2%, gi = 8% and Q = 40 is shown. Notice that these are the same width 
and parameters which allow trapping of a Townes soliton. 

We can see that with these parameters the stabilization is achieved although 
Eq. (23) does not predict trapping (in fact for these parameters p — 0). What 
is the difference with the stabilization of the Townes soliton shown in Fig. 
2?. A first important observation is that in the case of Gaussian initial data a 
readjustment is produced soon by ejecting a significant part of the wave packet 
far from the trapping region. This corresponds to the first 10 time units where 
the width significantly increases due to the contribution of the outgoing wave. 
When this wave hits the absorbing region it is dissipated leading to the step 
in the norm evolution shown in Fig. 5(e). What remains trapped is indeed a 
Townes soliton. 



Thus, the use of initial data different from a Townes soliton leads to a split- 
ting of the solution into the soliton itself, which is the structure which can be 
stabilized, plus a certain amount of radiation which goes far from the region 
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Fig. 5. Results of numerical simulations of Eq. (3) showing stabilization of Gaus- 
sian initial data u(x, 0) = (l/^/vrXo)e- r2 / 2X o with X = 1.09 for parameter val- 
ues go = —2tt, g\ = 8tt, = 40. (a-c) Plots of |it(x, t)\ 2 on the spatial region 
[—2, 2] x [—2, 2] corresponding to times of maximum (a,c) and minimum (b) width 
of the solution, (d) Evolution of the width X(t) and of the amplitude max \u(-, t)\. 
(e) Evolution of the norm r)(t). 

of interest. Because of this fact one must be careful to eliminate the outgoint 
part of the radiation in the numerical simulations. We have verified that when 
the absorbing region is absent the numerical simulations are misleading and 
the trapping effects are drastically altered (in fact, when zero boundary con- 
ditions at a given distance are used, there appear reflections in the boundary 
and essential changes on the results take place which would lead to spurious 
destabilization of the system after very short trapping times). 

It is known that Eq. (3) with a cubic nonlinearity and two spatial dimen- 



n(0i 

n no 
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sions corresponds to the so-called critical case [2,8] in which diffraction and 
self-focusing are nearly balanced and collapse is extremely sensitive to per- 
turbations and to changes in the initial conditions. For this reason although 
Gaussians look roughly like the Townes soliton they are not able to capture 
the delicate balance between diffraction and nonlinear focusing present in this 
case. When a Townes soliton is taken as a initial condition in Eq. (3) there is 
only a very small outgoing component which spreads out (due to numerical 
errors) and the system responds to the parametric perturbation as a whole. 
This fact implies that the moment equations may describe reasonably well the 
dynamics of the system only near the Townes soliton, since they only deal 
the global dynamics of the system. However, when the initial condition is a 
Gaussian, the component which spreads out cannot be captured neither by a 
moment-type formalism nor by the usual variational methods. 



5 Three-dimensional systems 



5. 1 Analysis of the moment equations 



Let us now consider the case n — 3, then Eq. (10) reads 

* = (24) 

Let us first assume that g is continuous and T-periodic, then ifg = 7pf T g>0 
bound states cannot exist. The reason is that, by Massera theorem, a bound 
state would imply the existence of a periodic solution. But if there is a periodic 
solution X, multiplying (24) by X 4 and integrating over a period we get 
> -4 / T X 3 X 2 - Qi / T X = gQ 2 T > 0, which is a contradiction. Therefore, 
a necessary condition for the existence of bounded solutions is 



9 rj-y 



fIo 9< °- (25) 



Note that this result is the opposite to the one inferred in Ref. [21] from direct 
numerical simulations an equation similar to (but restricted to the class of 
gaussian initial data) that of Eq. (24). 

Our goal is to find stable periodic solutions of Eq. (24) in three dimensions. A 
first useful observation is that an arbitrary positive T-periodic function X(t) 
is a solution of (24) if g(t) = Q^ 1 (X 4 X — QiX). In the following result, we 
use the definition of positive part of a function as f(t) + = max{0, f(t)}. 

Theorem 1 Let X(t) be a T-periodic positive function such that 
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(n) T f 
Jo 



T X + 

T <1 



Then, there exists g it) T -periodic function such thatX(t) is aT -periodic linear 
stable solution of (24). 

Proof. If we choose g(t) = Q% l (X 4 X — Q\X), then we know that X[t) is a 
T-periodic solution of (24). The linearized equation is 

Now, hypotheses {i) — (ii) correspond to the assumptions of the classical Lya- 
punov's criterion for linear stability (see for instance Theorem 1 in [41]). 

In many cases, hypotheses (i) — (ii) can be verified numerically. For instance 
taking g[t) = — Q^ 1 ^ + sin t) 4 sin £ — Q^ 1 (5i(15 + sint), then equation (24) 
has the solution X(t) = 15 + sint which is linearly stable. 

A different way to build solutions to Eqs. (24) is to choose g(t) = G(t)X(t) 
then Eq. (24) becomes equivalent to Eq. (11) and thus all the results of exis- 
tence of breathers in two dimensions can be applied. 

In conclusion, from the analysis of Eq. (24) we get that trapping could be 
possible in a three-dimensional scenario provided the dynamics of the system 
is well described by the moment equations (24). 



5.2 Failure of moment equations for symmetric 3D systems 



To check the validity of the predictions made on the basis of Eq. (24) we have 
compared its predictions with the behavior of the solution as obtained from 
direct numerical simulations of Eq. (3). We have used as initial data solutions 
of Eq. (19) for some \i and g s values obtained with a shooting method as 
we described above for the two-dimensional case. When n = 3, the relations 
between soliton solutions of different widths are 

R,{r) = ^R^r), X, = X 1 /^\ (26) 

The norm is not conserved by the previous change and verifies the relation 
•Hp = Vi/t* 1/A - 

In all the simulations we have taken g(t) — g + gi cos(fii) with parameters 
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go, gi and f2 chosen to give stable behavior on the basis of Eq. (24), 



Our numerical results show that the stabilization predicted by Eq. (24) does 
not occur for solutions of Eq. (3). Although it is possible by a suitable choice of 
go, gi and Q to change the behavior from collapsing to expanding or viceversa 
from the behavior which corresponds to the g\ = case, we were not able to 
find stable breather solutions for any choice of the parameters. In the region 
in which the solutions change from collapsing to expanding it is possible to 
fine-tune the parameters to get a solution which holds for some time, but this 
structure is unstable and finally the solution loses its profile. 

Thus, the moment equations fail to describe the dynamics of the system in 
three dimensions. The reason is simple: when the nonlinear term in Eq. (3) 
is cubic and n = 3 we are in the so-called supercritical case in which collapse 
occurs by formation of a localized spike on a finite-amplitude background. This 
behavior cannot be captured neither with the quadratic phase approximation 
[Eq. (7)] nor with any type of variational ansatz. 



5.3 Three dimensional systems with confinement along one direction 

Although stabilization of spherical structures seems not possible in three- 
dimensional scenarios, we could think about the situation in which the three- 
dimensional system has cylindrical symmetry. It would seem plausible that if 
the system is coin-shaped stabilization might take place because, in fact, this 
system is close to a quasi two-dimensional one and, therefore, the results in 
two dimensions could be applied. In fact, a numerical simulation reported in 
Ref. [20] supports this conjecture. Here we will try to make a more systematic 
analysis of the phenomenon. 

To get a qualitative understanding of the phenomenology for this situation let 
us consider Eq. (2). We will take gaussian functions as initial data and obtain 
the evolution equations for the parameters by the use of the collective coordi- 
nates method [12] to find the evolution of the widths in each spatial direction. 
The idea is to restate the problem of solving Eq. (2) as a variational problem, 
corresponding to a stationary point of the action related to the Lagrangian 
density C. So, the problem is transformed into the problem of finding i/j(x.,t) 
such that the action 



is extreme. This problem is as complicated as solving the original NLSE. The 
idea of the method is to restrict the analysis to a set of trial functions. One 




(27) 
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possible choice is to take a Gaussian ansatz of the form 

1 n 

J.l-X f\ - TT -x 2 j /2w]+ix j a j (t)+ixp j (t) 

^/^n =iwAt)V2 Yie , (28) 

which verifies the following relations 



Wj 



(t) = ^ \u\ 2 d n x^j = 1, (29a) 
(t) = 2 y \u\ 2 x 2 j( Tx. (29b) 



The variational method leads to the following evolution equations for the 
widths 

In the three-dimensional case the equations are 



a2/ , 1 #(t)/(27r) 3 / 2 , x 

U>1 w(w 2 W 3 

a2/x 1 #(t)/(27r) 3 / 2 . oiM , 

ui 2 + A 2 (t )w 2 = ^ + \ , 31b 

.. x2/ . 1 (?(t)/(27T) 3 / 2 

^3 + ^M = ^+ 2 • 31C 

If we consider cylindrically symmetry solutions we have u>i = w 2 = w, this 
width w being equivalent to the width X arising in the radially symmetric 
two-dimensional case. 

X= (/ HV + ^) 1/2 =(^) 1/2 =». (32) 

Therefore, if we consider the case where the potential is present only along the 
z direction and we suppose cylindrical symmetry the equations for the widths 
are 



.. 1 ff(t)/(27T) 3 / 2 

w = —-\ - , (33a) 

(33,, 

Our goal is to stabilize some solutions of this model, i.e. to get periodic solu- 



19 



tions for (w(t),W3(t)). In particular, if we impose that 

9(t) = 92p(t) 
(2ir) 3 / 2 w 3 (t) 2vr ' 



(34) 



where (?2d(£) is an specific modulation of the nonlinear term allowing trapping 
in two dimensions, then Eq. (33a) is the same equation that in the radially 
symmetric two-dimensional case (Eq. (23)) and, therefore, by taking the three- 
dimensional modulation of the form of Eq. (34) we will have stabilization of 
the width w. Now, Eq. (34) can be written as 

g{t) = (27T) 1 / 2 g 2D (t)w 3 (t), (35) 

and the problem arises when we realize that w 3 (t) is not known, but it evolves 
according to Eq. (33b), so Eq. (35) is useless in practical cases. Nevertheless, 
we could think about the possibility of stabilizing w 3 around the equilibrium 
point of Eq. (33b) if the nonlinear term were absent. This equilibrium point 

— 1/2 

is Wzfi = A 3 .To ensure that the value of u> 3 is as similar to the equilibrium 
value as possible, we have to impose that 



for all t, and this implies that 

92D(t) 



« -T> (36) 

w 3,0 



A 3 max 



2nw 2 (t) 



(37) 



If Eq. (37) were satisfied we could take w 3 (t) ~ w 3i o = A 3 in Eq. (35) 
and from this equation we get a prediction for the values of the modulation 
we should take to get stabilization. As we know the values of g^D and w to 
stabilize in two dimensions, Eq. (37) is the condition to get stabilization in 
three dimensions (with the potential along z). As we said before, this condition 
means that the more coin-shaped the system is, the better the stabilization will 
be, leading to a quasi two-dimensional system. Again, the variational method 
predicts that stabilization can occur only if g$ < (2tt) 1 ^ 2 (—2-k) but we will see 
later that it is also possible for other g values. 

In Fig. 6 we plot the evolution of the widths after solving numerically Eqs. (33) 
for different values of A 3 . We have taken g = (2n) 1/2 (-2.2n), g 1 = (2n) 1/2 8n, 
Q = 40, w(0) = 1.09 and w 3 (0) = w 3t0 = A 3 . We expect that the values of 
w during the evolution will be about w — 1, as in 2D case, so the condition 
(37) says that A 3 3> 5, approximately. We can see that the greater the value 
of A 3 is the better the stabilization is, according to our previous estimates. 

We have compared these results with the simulations of the full NLSE for the 
case when the potential is restricted to the z-axis and the system is strongly 
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Fig. 6. Stabilization of a Gaussian according to numerical simulations of Eq. (33) 
with parameters g = (27r) 1 / 2 (-2.27r), gi = (2vr) 1 / 2 87r, Q = 40, w(0) = 1.09, 
w 3 (0) = W3fi = A3 1/2 

is shown for (a) A 3 = 60, (b) A 3 = 80 and (c) A 3 = 100. 



3 . The evolution of the widths w (dashed) and W3 (solid) 



coin-shaped. To do so we have used a pseudo-spectral fully asymmetric 3D 
evolution method combined with a second order split-step method to advance 
in time. Typical grid sizes were of 128 x 128 x 64. 

In Fig. 7 we present the results for g = (27r) 1/2 (-27r), g x = (2n) 1/2 8n, Q = 
40, A 3 = 100 and as initial data we take a Gaussian with w(0) = 1.09 and 

— 1/2 

^3(0) = A 3 . We see that, as in the two-dimensional case, part of the wave 
packet goes outside and after a readjustment the solution oscillates in the same 
way that in two dimensions. Moreover, the width w 3 remains nearly constant 
during the evolution. 



6 Conclusions and discussion 



In this paper we have provided a deeper understanding of the phenomenon 
of stabilization of solitons of the nonlinear Schrodinger equation by means of 
the control of the nonlinear term. We have developed the moment equations 
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Fig. 7. Results of numerical simulations of Eq. (2) showing stabilization of Gaus 
sian initial data ^(x,0) = (l/vrV^O^O) 1 / 2 ^ 



A. 



-1/2 



-(x 2 +y' 2 )/2w(0) 2 e - 



for parameter values go 



z 2 /2w 3 (0) 2 with 

(2^) 1 /2 ( _2vr), 



widths w(0) = 1.09 and w 3 (0) - 
gi = (27r) 1 / 2 87r, 9, = 40 and A 3 = 100. (a-c) Isosurface plots of \u(x,t)\ 2 for 
\u\ 2 = 0.01 corresponding to times of maximum (a,c) and minimum (b) width w(t) 
of the solution, (d) Evolution of the widths w(t) and ws(t) and of the amplitude 
max|n(-,t)| of the solution, (e) Evolution of the norm ij(t). 

for 2D and 3D systems and obtained, on the basis of their rigorous analysis, 
precise conditions for the stabilization of 2D systems. 



Taking as initial data Townes solitons or Gaussians for simulations based on 
Eq. (3) we have shown that the former is the structure which is stabilized and 
that other initial data which can be trapped must eject a fraction of the wave 
packet in the form of radiation to accomodate to this specific solution. 

Also we have analyzed the three-dimensional situation. Here we have made an 
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extensive search of stable regions according to our moment equations improv- 
ing and extending the analysis of Ref. [21]. We have justified why moment-type 
equations cannot be used to predict the dynamics of the system. Only limited 
time stabilization is possible when the parameters are fine-tuned to very pre- 
cise values. Finally, when a strong trapping along one specific direction is kept 
the system becomes effectively two-dimensional and it can be described again 
by variational methods whose predictions agree well with the full numerical 
simulations of the problem. In the latter case three-dimensional confinement 
is possible although there are only two spatial directions along which the solu- 
tion is trapped by the nonlinear forces plus the stabilization mechanisms while 
the other direction is trapped by harmonic forces provided by the potential. 

It is remarkable that the nonlinear Schrodinger equations support these sta- 
bilized structures which open new fields for applications, in fact these are 
the first stable structures obtained in the framework of the cubic nonlinear 
Schrodinger equation. 

Many extensions of this work are possible. First, it would be very interest- 
ing to study the robustness of these breather-type solutions under different 
perturbations, e.g. by mutual collision of different structures in single and 
multicomponent systems. Secondly it could be interesting to try to stabilize 
other stationary structures different from the two-dimensional Townes soliton. 
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